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Abstract 

We examine the Galilean invariance of standard lattice Boltzmann methods for 
two-phase fluids. We show that the known Galilean invariant term that is cubic in 
the velocities, and is usually neglected, is the main source of Galilean invariance 
violations. We show that incorporating a correction term can improve the Galilean 
invariance of the method by up to an order of magnitude. Surprisingly incorporating 
this correction term can also noticeably increase the range of stability for multi- 
phase algorithms. We found that this is true for methods in which the non-ideality 
is incorporated by a forcing term as well as methods in which non-ideality is directly 
incorporated in a non-ideal pressure tensor. 
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1 INTRODUCTION 



The issue of Galilean invariance for lattice Boltzmann methods, and in partic- 
ular for non-ideal fluid simulations based on an input pressure tensor [1], has 
received considerable attention in the past [2,3,4]. It was noticed that these 
methods had a severe problems with Galilean invariance and careful expan- 
sion methods elucidated a set of correction terms for the pressure tensor to 
improve Galilean invariance. The attraction of methods based on an input 
pressure tensor is that this pressure tensor can be easily derived from an in- 
put free energy. This immediately delivers predictions for the phase-diagram, 
the interface profiles and the surface tension. 

Others held that the problem lay deeper, and that it was inappropriate to in- 
clude the non-ideal terms for the method in an input pressure tensor, and that 
it is appropriate to put the non-ideal terms in a Vlasov-like forcing term[5,6]. 
Unfortunately, it has been difficult to relate the approach by Shan and Chen [5] 
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to a free energy. Meanwhile others have taken to use a forcing term that cor- 
responds to the divergence of the pressure tensor [6]. This should, in principle, 
lead to an equivalent model to the pressure tensor approach. We will show 
later, that there are subtle differences though. On the up-side, using a forc- 
ing approach does not require similar correction terms to the pressure tensor 
approach, and therefore it has been labeled "Galilean invariant" . 

Recently we have been working on including Lees-Edwards-like boundary con- 
ditions in lattice Boltzmann [9] to simulate sheared systems. It turns out that 
these boundary conditions are very sensitive to Galilean invariance violations. 
None of the above methods seemed to show a sufficient level of Galilean in- 
variance which lead us to this closer investigation of the origins of Galilean 
invariance in lattice Boltzmann methods. 

In a lattice based method the lattice represents a fixed reference frame, and it 
is not surprising that this should show up to some order. So we measured the 
error for a variety of situations and tried to quantify the main contribution. 
We found that the well known V 2 (pw 3 ) term in the Taylor expansion of the 
momentum equations, which is usually neglected, is responsible for the major- 
ity of the error. We then demonstrate how including an additional correction 
term can significantly increase the Galilean invariance of lattice Boltzmann 
methods. 



2 The lattice Boltzmann method 

We can write the lattice Boltzmann method in a general way that neatly 
separates the ideal gas contributions from the non-ideal contributions. To do 
this we write 

fi(r + vAt, t + At)- fi(T, t) = -(/°(x, t) - /i(x, t) + d) + F t (1) 

T 

where ff is the contribution function for an ideal gas, the Fi are the contribu- 
tions of a forcing term and the G{ allow us to manipulate the pressure tensor. 
As usual the moments for the ideal-gas distribution function are 

fi = Pi fi V ia =PUai fiViaVip = (? s pb afi + pU a Up, 

i i 

i 6 

where p = J2i fi is the mass density, u = J2i fi v i/p is the mean fluid velocity 
before the action of the forcing term Fj. The velocity of sound is given by c s = 



2 



l/y/3. For all models with vf a = 1 (like D1Q3, D2Q7, D2Q9, D3Q15, D3Q19 
or D3Q27) we have vf a = v ia which means that e.g. Y,i fi v f x = Hi fi v ix = pu x 
making a correction term Q necessary for the third moment. One usually 
chooses Qafo = —pUaUpu^. It is this term that will lead to the leading Galilean 
invariance problems. 

To simulate fluids with a non-ideal equation of state we can introduce either 
a forcing term Fj with the moments 

Yl F i = °> F i V ic = P a ai F i V iaVi(S = p(a a llp + ClpU a ), 

i 

1 

FiV ia v ip v h = -{pa a 5 Pl + pa p 5 a p + pa 7 5 a/3 ) 
or a pressure term Gi with the moments 

J2Gi = 0, G i V i* = °> Yl G i V ioVifi = A a p, Gv ia Vii3V h = 0. 

i i 

Here pa is a forcing term and A is a pressure term. We will see below how 
these terms can be used to introduce non-ideal terms into the equation of state 
or simply to correct the above mentioned deficiencies of the velocity set. 

The knowledge of these moments is sufficient to perform a Taylor expansion 
(or equivalently a Chapman Enskog multi-scale expansion) of equation (1). If 
we define the macroscopic velocity u as u = u + a/2 we obtain the continuity 
equation 

d t p + d a (pu a ) = (2) 

and a momentum conservation equation 

d t {pu a ) + dp{pu a up) = -d/3(pc 2 s 5 a/ 3 + A af3 ) + pa a 

+dp [vp{dpu a + d a up + <9 7 w 7 o" a/3 )] (3) 
—vdpluadyAfa + upd~fA ai + d p A al3 d^(pu^) + <9 7 Q a/ 3 7 ] 

where the kinematic viscosity is v = (r — 1/2) c^. This equation is the Navier- 
Stokes equation, except for the terms in the last line. The condition of Galilean 
invariance requires that a description in a reference frame S another reference 
frame S' translating with a constant velocity uq be related by a spatial coor- 
dinate transformation x = x' — u^t and a translation of velocities u = u' — uq. 
It is easy to see that only the terms in the last line of eqn. (3) are not Galilean 
invariant. 
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Let us first consider (3) for an ideal gas. In this case a = and A = so the 
only non-Galilean invariant term to second order is Q. We can eliminate this 
error term by introducing a well crafted forcing term of the form 



We should mention here that we can avoid this problem and choose Q = 
if we use a velocity set that is large enough [7]. But for now we want to stay 
with the standard velocity sets. 

To simulate a non-ideal system we want to obtain a pressure term in the first 
line of (3) that is the divergence of a pressure tensor derived from Thermody- 
namics. In particular we need 



to first order. For the forcing approach we choose A = and this equation 
defines a. We will refer to this approach as "Forcing" in our comparisons. To 
improve Galilean invariance for this approach we can add the same additional 
forcing term (4) to restore Galilean invariance to second order. This approach 
we will refer to as "ForcingQ" . 

If we chose to introduce the non-ideal term in A we can recover the algorithm 
of Swift et al.[l]. This corresponds to choosing A a p = P a p — pc 2 s and a = 0. 
We refer to this approach in the following as "Pressure" . 

This approach is deficient in as much as the correction terms in the third 
line in (3) are unphysical and will lead to severely non-Galilean invariant 
behavior. This problem was addressed by Holdych et al. [2] and independently 
by Inamuro et al. [3] . Because the divergence of the pressure tensor is small to 
first order, the main contribution to the error terms comes from the density 
alone. Therefore choosing 

A al 3 = PafS ~ pC 2 s p5af3 ~ v(d a pUp + d p pU a + d^pU^) (6) 

will only leave error terms of the order d 2 P which can be assumed to be 
small in systems close to equilibrium. This approach will be referred to as 
"Holdych" . The restriction of being close to equilibrium was late revisited by 
Kalarakis et al.[4] who suggested an improvement to these corrections. 

However, the Holdych approach as well as all other later approaches are still 
deficient in that they do not correct the Q term. Clearly there is a multitude 
of possible combinations of choices for A and a that will lead to a Galilean 
invariant form of eqn. (3). One other choice we examined is 



pa a = vdpdjQafa 



(4) 




(5) 
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A a p = P af3 - pc s , 

a a = -dpivd^u^Pfo - pc 2 s 5^) + up(A ai - pc 2 s 5 ai ) + (d p A afS - c 2 s SaP)d J (pu J )]}. 

which has in common with the Kalarakis[4] approach that it is not limited to 
systems close to equilibrium. The above choice leads to a Galilean invariant 
momentum equation given by 

d t (pu a ) + d p (pu a up) = -dp(pc 2 s 5 a p + A a/3 ) + pa a 

+d p [u(P^d 7 u a + Pa-yd^up + P a/3 <9 7 w 7 ))] (7) 

This scheme has a tensorial interface viscosity. We refer to the resulting 
method as "PressureQ" . 



3 The non-ideal gas 

For simplicity we will consider a non-ideal gas with a 4 -free energy [8]. For 
such a system we can calculate the phase diagram and the surface tension 
analytically, simplifying the analysis. 

For a critical density p c , critical temperature T c and critical pressure p c we 
obtain for the pressure tensor 

Pa/3 = [Pc(0 + 1) W - 20 + 1 + 29) - npV 2 p - ^{Vp) 2 ]5 aP + Kd aP dpp{%) 

where = (p — p c )/p c is the reduced density and 9 = (3(T — T c )/T c is the 
reduced temperature, where (5 is an arbitrary constant. The equilibrium values 
for the density are given by 

p° = p c ± V^9. (9) 

The definition of the pressure tensor is all that is needed to define the lattice 
Boltzmann methods for non-ideal fluids as we explained in section 2. 

We performed simulations of an equilibrium system containing one domain of 
gas and one of liquid at different imposed velocities. As mentioned above, while 
the different approaches are very similar as far as the expansion to second order 
is concerned, there are noticeable differences in the behavior of the methods. 
Firstly let us compare the numerical results for the phase-diagram shown in 
Figure 1. On the one hand we notice that the ability of the pressure based 
methods to reproduce the analytical phase-diagram is noticeably better than 
the corresponding forcing method. On the other hand we we see that the range 
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Fig. 1. Phase diagram for p c =0.1, k = 0.1, n c = l,p c = 0.1 and v = 1/6 for the 
"ForcingQ" and "Holdych" approaches. Note that the Holdych approach is better 
at reproducing the analytical phase diagram where as the ForcingQ approach has a 
larger range of stability. 
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Fig. 2. Deviation from the analytical density difference for the different methods as 
a function of uO. k = 0.1, n c = 1.0, p c = 0.42, (3 = 0.1, 6 = -0.03 and v = 1/6. The 
data end at values of uq at which the methods became unstable. 

of stability for the forcing method is larger, leading to a larger maximum 
density ratio (of about six) which can be simulated with this method. 

Next we examine the effect of an imposed velocity on the equilibrium den- 
sities. In Figure 2 we show the deviation of the predicted liquid-gas density 
difference from the measured one. We notice that the original pressure method 
has excessive deviations even for small u D . We were surprised to see that in 
these simulations the gas and liquid domains were actually stationary, even 
with an imposed velocity uq. The domains made up for the imposed veloc- 
ity by evaporation and condensation mechanisms and a faster velocity in the 
gas than in the fluid. We see that the error is generally less for the corrected 
pressure methods than for the forcing methods. 

We were very surprised to see that the Q corrections did not only improve 
that Galilean invariance of the method, they substantially increased the range 
of stability. This is true both for the range of stable velocities uo as well as 
the accessible density ratios (even at u = 0. This was an unexpected benefit 
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Fig. 3. Ed(uo) for simulation parameters k = 0.01, n c = 1, p c = 0.42, 6 = —0.03, 
iteration= 10000 for the different methods. 

of our study and we still do not understand why the correction term has such 
beneficial effects on the stability. 

We now need to quantify the Galilean invariance error for the advection of an 
interface profile. For the velocity the analytical solution is a constant, so we 
can define an error function 



fiflW^EW 1 )-*' (io) 



This measure is effectively time independent, as is appropriate for this equi- 
librium consideration. 

In Figure 3 we show the Galilean invariance error E^{u ) for the different 
methods we described. As mentioned above the original pressure approach 
performs very poorly, in fact refusing to advect the domains relative to the 
lattice. The Holdych approach improves on this significantly leading to an 
advection of the profile. We were surprised by the non-monotonic behavior 
of the error, leading to a minimum at u$ = 0.2. We will come back to this 
later. Adding the Q correction to the Holdych approach leads to a noticeable 
improvement for velocities as small as Uq = 10 -3 . The Forcing approach leads 
to a good behavior at small u but increases rapidly with u . Its behavior is 
significantly improved for uq > 0.03 by including the Q correction in the in 
the method. 

When interpreting the above results, it is important to remember that the 
parameter space for the Galilean invariance problem includes not only uq 
but also the parameters determining the equilibrium density profile k, 9 and 
p c as well as the relaxation time r. This parameter space is so large as to 
make it nearly impossible to examine it exhaustively. But we want to discuss 
at least the dependence on k, which is related to the interface width and 
the surface tension. It is also important to look at this when one wants to 
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Fig. 4. Ed(uq) as a function of n for different values for uq for the PressureQ method. 
Other parameters are n c = 1, p c = 0.42, = —0.03. 

fairly compare the pressure an forcing approaches. It turns out that while 
the pressure approach reproduces the analytical interface profile fairly, the 
forcing approach does not, at least not for the nominal value of k. The forcing 
approach leads to a much wider interface so that it would be fairer to compare 
the forcing approach to a pressure approach with a larger k. In Figure 4 we 
can see that this matters a lot. The Galilean invariance error Ed is largest for 
small values of k corresponding to thin interfaces. And for small values of k 
this absolute error decreases much slower with decreasing uq. This is probably 
the reason that the pressure approach performed much worse than the forcing 
approach (with the same nominal k value) for small uq. 

A close examination of Fig 4 also shows that the non-monotonic behavior of 
E d {uq) is related to the small k behavior. The graph suggests that the error 
function becomes monotonic for large k. 



4 Summary 

We have shown that the usually neglected error term in the Navier Stokes level 
momentum equation derived for standard lattice Boltzmann methods leads to 
noticeable Galilean invariance violations. These violations are noticeable even 
in the range of usual velocities of u < 0.1, but become dominant for larger 
velocities. 

A carefully defined forcing term can remove the non- Galilean invariant terms 
recovering the Navier Stokes equation. We have shown that this approach is 
also effective in practice, reducing the Galilean invariance error substantially. 
The correction term also had the added benefit of increasing the range of 
stability for the multi-phase applications, leading to a larger range of stable 
velocities and even, perhaps surprisingly, to a larger range of density ratios 
that can be simulated by the methods. 




i 



8 



References 



[1] M.S.Swift,W.R.Osborn, J.M. Yeomans,P%s. Rev. Lett. 75, 830 (1995). 

[2] D.J. Holdych,D. Rovas,J.G. Georgiadis, and R.O. Buchius, Int. J. of Mod. Phys. 
C 9, 1393 (1998). 

[3] T. Inamuro, N. Konishi, F.Ogino, Corny. Phys. Comm. 129, 32 (2000). 

[4] A.N. Kalarakis, V.N. Burganos, and A.C. Payatakes, Phys. Rev. E 65, 056702 
(2002). 

[5] X. Shan,H. Chen,Phys. Rev. £49, 2941(1994). 

[6] L.-S. Luo, Phys. Rev. E 62, 4982 (2000), L.-S. Luo, Phys. Rev. E 62, 4982 
(2002). 

[7] Y.-H. Qian, Y. Zhou, ICASE Report No.98-38, NASA/CR-1998-208701(1998). 
[8] A.J. Briant, A.J. Wagner, and J.M. Yeomans, Phys. Rev. E 69, 031602 (2004). 
[9] A.J. Wagner and I. Pagonabarraga, J. Stat. Phys. 107, 521 (2002). 



9 



